retaildata <- read.csv("retail.csv")
mytimeseries <- ts(retaildata[,4], frequency=12, start=c(1982,4))
autoplot(mytimeseries)There is a strong trend to around 2008, weak seasonality (with a jump in December each year), and some evidence of cycles (with dips in 2001 and 2012).
bicoalbicoal is annual data, so you can’t do the seasonal plots.
doleOther series are similar.
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 32.269, df = 8, p-value = 8.336e-05
##
## Model df: 0. Total lags used: 8
There is some remaining autocorrelation in the residuals: the Null of no joint autocorrelation is clearly rejected. We can also see a significant spike on the seasonal (4th lag) in the ACF. There is considerable information remaining in the residuals which has not been captured with the seasonal naïve method. The residuals do not appear to be too far from Normally distributed.
train <- window(mytimeseries, end=c(2010,12))
test <- window(mytimeseries, start=2011)
autoplot(cbind(Training=train,Test=test))## ME RMSE MAE MPE MAPE MASE
## Training set 5.502402 27.74935 19.33784 3.369450 10.447161 1.0000000
## Test set 5.227586 23.40458 18.60000 1.619239 8.172045 0.9618449
## ACF1 Theil's U
## Training set 0.8703252 NA
## Test set 0.4544630 0.9518744
The number to look at here is the test set RMSE. That provides a benchmark for comparison when we try other models.
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 1256.8, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
The residuals do not look like white noise there are lots of dynamics left over that need to be explored. They also do not look close to normal, with very long tails.
The accuracy measure are always sensitive to the training/test split. There are better ways to check the robustness of the methods in terms of accuracy such as using a tsCV().
e_mean <- tsCV(mytimeseries, meanf)
e_naive <- tsCV(mytimeseries, naive)
e_snaive <- tsCV(mytimeseries, snaive)
e_drift <- tsCV(mytimeseries, rwf, drift=TRUE)
# Construct squared CV errors matrix
e2 <- cbind(Mean=e_mean^2, Naive=e_naive^2,
SNaive=e_snaive^2, Drift=e_drift^2)
# Remove rows with any missing for a fair comparison
e2 <- na.omit(e2)
# Find MSE values
colMeans(e2)## Mean Naive SNaive Drift
## 4685.8665 410.0438 691.7144 411.4290
## ME RMSE MAE MPE MAPE MASE
## Training set -2.740622 26.56395 19.38206 -2.872958 10.05648 0.9560879
## ACF1
## Training set -0.005983602
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
## ACF1
## Training set 0.01348202
## ME RMSE MAE MPE MAPE MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
## ACF1
## Training set -0.003195358
These RMSE values are not really comparable because the models have different numbers of parameters.
The best model is the last one.
##
## Ljung-Box test
##
## data: Residuals from Damped Holt's method
## Q* = 10.782, df = 5, p-value = 0.05587
##
## Model df: 5. Total lags used: 10
The residuals are pretty good and (just) pass the Ljung-Box test. However, the forecast make little sense (with the prediction intervals going negative very quickly).
The seasonal variation increases with the level of the series. So we need to use multiplicative seasonality.
## ME RMSE MAE MPE MAPE MASE
## Training set -0.2054187 9.280748 6.244629 -0.2790438 3.338901 0.3336358
## ACF1
## Training set 0.05182697
## ME RMSE MAE MPE MAPE MASE
## Training set 0.5050393 9.194427 6.128123 0.2458937 3.290157 0.3274111
## ACF1
## Training set -0.04640084
There is not much difference between these models, and we would expect the trend to continue, so I would prefer to use the non-damped version.
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 41.483, df = 8, p-value = 1.693e-06
##
## Model df: 16. Total lags used: 24
There are significant correlations in the residuals, especially at lag 12. So these residuals do not look like white noise.
train %>% window(end=c(2010,12)) %>%
hw(seasonal='multiplicative', damped=FALSE) %>%
accuracy(x=mytimeseries)## ME RMSE MAE MPE MAPE MASE
## Training set -0.1212843 9.260807 6.054543 -0.03697659 3.442954 0.3130931
## Test set 10.1364618 19.695875 15.215508 4.30118709 6.977221 0.7868257
## ACF1 Theil's U
## Training set 0.09051818 NA
## Test set 0.61161198 0.909366
The test set RMSE is much better than for the seasonal naïve method.
## ETS(M,Ad,M)
##
## Call:
## ets(y = mytimeseries)
##
## Smoothing parameters:
## alpha = 0.7718
## beta = 0.0025
## gamma = 0.0561
## phi = 0.98
##
## Initial states:
## l = 64.0899
## b = 0.2813
## s = 0.9869 0.9362 1.0043 1.1705 1.0216 1.0293
## 0.9848 0.9998 0.993 0.9436 0.9694 0.9606
##
## sigma: 0.0464
##
## AIC AICc BIC
## 4404.041 4405.697 4477.272
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.5339249 9.613425 6.241919 0.1970991 3.257992 0.333491
## ACF1
## Training set -0.03180379
This is equivalent to a damped-trend multiplicative Holt-Winters’ method. The very small \(\beta\) and \(\gamma\) values show that the seasonality and trend change slowly. There is also very little damping (\(\phi\) is close to 1).
Not too bad given the noisy data.
These are great forecasts which have adapted to the trend and seasonality really well.
The seasonality seems to be over-done, but the wide prediction intervals are probably sensible apart from the fact that they go negative.
These data are cyclic, and ETS does not handle cycles in time series.